/*
 heat-spectral lib. solves time-dependent PDEs using SEM-NI method.
 Copyright (C) 2009, Altieri Michele and Stefanoni Andrea
 Copyright (C) 2008, Proverbio Alessandro and Radice David

 This program is free software: you can redistribute it and/or modify
 it under the terms of the GNU General Public License as published by
 the Free Software Foundation, either version 3 of the License, or
 (at your option) any later version.
 
 This program is distributed in the hope that it will be useful,
 but WITHOUT ANY WARRANTY; without even the implied warranty of
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 GNU General Public License for more details.
 
 You should have received a copy of the GNU General Public License
 along with this program. If not, see <http://www.gnu.org/licenses/>.
*/

/**
 * \file problem.hpp
 * \brief This file contains the Problem class
 *
 * This class is to be used for solving the parabolic problem with a SEM-NI
 * method. It is designed for rectangular domains and for any kind of
 * boundaries conditions.
 */

#ifndef PROBLEM_HPP
#define PROBLEM_HPP 1

#include <cassert>
#include <fstream>
#include <list>
#include <valarray>
#include <vector>

#include "mesh.hpp"
#include "point.hpp"
#include "typedefs.hpp"
#include "umfmatrix.hpp"

namespace spectral {

/**
 * \class Problem
 * \brief Assemble SEM-NI parabolic problem
 *
 * The class configure a parabolic problem over a rectangular domain with
 * Dirichlet or Neumann boundaries conditions and with any kind of
 * bilinear form and functional, even time-dependent.
 * \param BIL Bilinear form, the BIL class must be an object like this
 * \code
template<class MESH>
class Bilinear1 {
    public:
        Bilinear1(MESH const & me);

        template<class U, class V> inline
        double operator()(U const & u, V const & v, int n) const;

        // Boundary conditions type
        bctype bcond[4];

        // Dirichlet conditions (as a function of the point)
        struct dirichlet{
            double operator()(double const & x, double const & y) const;
        };
 
        // Neumann conditions (as a function of the point)
        struct neumann{
            double operator()(double const & x, double const & y) const;
        };
 
        // Methods on (which apply on different borders)
        template<unsigned int N, typename T=void>
        struct on : public dirichlet {};       // If a dirichlet condition is the
                                               // default type
 
        template<typename T>
        struct on<1,T> : public neumann ();    // If, for example, we want to impose
                                               // a neumann condition on border 1
    private:
        MESH const & me;
};
 * \endcode
 * \param F Functional, the F class must be an object like this
 * \code
template<class MESH, class FUNC>
class Functional {
    public:
        Functional(MESH const & me);

        template<class V> inline
        double operator()(V const & v, int n) const;
    private:
        MESH const & me;
        FUNC f;
};
 * \endcode
 */

/**
 * \fn Problem::Problem(Point<double> const & AP,
 *                      Point<double> const & BP,
 *                      int nx,
 *                      int ny,
 *                      int nnodes,
 *			double deltat,
 *			double theta,
 *                      int nmax,
 *                      double legDelta,
 *                      double legEps,
 *                      int order=1)
 * \brief The Problem class constructor
   \param AP Bottom-Left point of the dominium
   \param BP Top-Right point of the dominium
   \param nx Number of elements on the x-axis
   \param ny Number of elements on the y-axis
   \param nnodes Number of nodes for the quadrature in every element
   \param deltat Step for the time-advancing scheme
   \param theta Relaxation parameter of the time-advancing scheme
   \param nmax Maximum number of iteration for the Newton-Horner method
   \param legDelta Minimum increment in the newton formula used to compute the GLL nodes
   \param legEps Critical resiudual in the newton formula used to compute the GLL nodes
   \param order Is the order of the problem in the week form. For example it's 1 for the standard advection-diffusion problem.
*/

/**
 * \fn void Problem::initialize(UMFMatrix & A,
 * 				UMFMatrix & M,
 * 				UMFMatrix & K,
 * 				std::valarray<double> & b,
 * 				std::valarray<double> & bold,
 * 				std::valarray<double> & u,
 * 				std::valarray<double> & uold,
 * 				std::valarray<double> & residual,
 * 				std::valarray<double> & delta,
 * 				std::valarray<double> const & bc)
 * \brief Initializes the problem
 *
 * This will generate all the matrices and the array necessary to solve
 * the parabolic problem. Actually, the structures that are going to be
 * modified during the resolution of the problem are now initialized only.
 *
 * \param A Stiffness matrix of the associated elliptic problem
 * \param M Mass matrix of the parabolic problem
 * \param K Matrix related to the resolution of the time-advancing step
 * \param b Known term related to the functional and the neumann boundary terms
 * \param bold Array containing the known term at the previous temporal step
 * \param u Array that contains the solution at the actual time
 * \param uold Array that contains the solution at the previous temporal step
 * \param residual Array that contains the update of the solution for each step
 * \param bc Vector that contains the solution on the boundary edges
 */

/**
 * \fn void Problem::initialize()
 * \brief Initializes the problem
 *
 * This method simply calls the version with several input parameters.
 * The interesting fact is that the final user can call directly this
 * function from the main file, in order to avoid viewing all the
 * algebraic structures.
 */

/**
 * \fn void Problem::matrixGen(UMFMatrix & A) const
 * \brief Stiffness matrix generator
 *
 * This will fill the element of the matrix in UMFPack standard with
 * the values related with the problem data. It does not computes the
 * effect of the border edges, which will be evaluated with the
 * patchMatrix method.
 */

/**
 * \fn void Problem::massMatrixGen(UMFMatrix & M) const
 * \brief Mass matrix generator
 *
 * This will generate the mass matrix for the parabolic problem in
 * UMFPack standard. This matrix is simply obtained from the
 * discretization of the integral operator
 * \f[
 * 	\int_\Omega \psi_i \psi_j \,\mathrm{d}\Omega
 * \f]
 * where \f$ \psi_i,\psi_j \f$ are the Lagrange basis functions.
 */

/**
 * \fn void Problem::initGen(std::valarray<double> & u) const
 * \brief Initial data generator
 *
 * This will interpolate the analytic initial data given from the user
 * and will put the discrete result into the array \c u. In this way
 * the time-advancing scheme will start from this array.
 */

/**
 * \fn void Problem::forceGen(std::valarray<double> & b) const
 * \brief Known term generator
 *
 * This will fill the element of the array representing the known term
 * of the system that will be solved. This method will not insert the
 * terms related with the neumann boundary conditions.
 */

/**
 * \fn void Problem::forceGen(std::valarray<double> & b, double const & t) const
 * \brief Known term generator when the forcing term is time-dependent
 *
 * This will fill, at every temporal step, the element of the array
 * representing the know term of the system that will be solved.
 * This method is different from the stationary one, because it has to
 * be called into the solver with a different parameter \c t at each
 * time-step.
 */

/**
 * \fn void Problem::patchMatrix(UMFMatrix & A) const
 * \brief Matrix adapting due to boundary conditions
 *
 * This function will modify the A matrix in order to fit the dirichlet
 * boundary conditions where necessary. Actually neumann boundary
 * conditions will not modify the stiffness matrix but only the known term.
 */

/**
 * \fn void Problem::bcGen(std::valarray<double> const & bc,
			   std::valarray<double> & b) const
 * \brief Boundaries conditions generator
 *
 * This will set the boundaries conditions both in the vector of the
 * known term. It requires the value of the boundary condition, and the
 * type, if dirichlet or neumann.
 */

/**
 * \fn void Problem::edgeBc<unsigned int border>(double const * const & newBcEdge,
						 std::valarray<double> & b) const
 * \brief Single border boundary conditions generator
 *
 * This will set the boundaries conditions related to the single border passed
 * in the template form. This method is called by bcGen which replicates this
 * procedure on the several border edges.
 */

/**
 * \fn void Problem::kappaMatrixGen(UMFMatrix & A,
 * 				    UMFMatrix & M,
 * 				    UMFMatrix & K)
 * \brief Building and inversion of K matrix
 *
 * This will computes the element of the K matrix, that is the linear operator
 * that needs to be inverted in order to solve the problem. Moreover, this method
 * provides the LU factorization of the matrix in UMFPack form. In this way the
 * resolution of each temporal step is definitely more direct, beacuse the great
 * part of the computation has already been done.
 */

/**
 * \fn void Problem::solve()
 * \brief Solves a singular step of the parabolic problem
 *
 * This simply computes the problem's residuals multiplying by matrix A the
 * solution at the previous time step. Then, it multiplies the inverse K matrix
 * by the residuals in order to obtain the difference from the new solution.
 * In the end, it sums this difference to the old solution, and ultimates the
 * temporal step. It returns
 * \f[
 * 	\mathbf{u}^{k+1} = \mathbf{u}^k + \mathrm{K}^{-1}\mathbf{r}^k
 * \f]
 * where
 * \f[
 * 	\mathbf{r}^k = \mathbf{f}^k - \mathrm{A}\mathbf{u}^k
 * \f]
 */

/**
 * \fn void Problem::solve(double const & t)
 * \brief Solves a singular step of the parabolic problem
 *
 * This is the same method of \p solve, but it solves parabolic problems with a
 * time dependent forcing term. In this case the function takes a \p double as
 * the time corresponding to the step \f$ k+1 \f$ of the iterative procedure.
 * So the method returns
 * \f[
 * 	\mathbf{u}^{k+1} = \mathbf{u}^k + \mathrm{K}^{-1}\left(\mathbf{r}^k
 * 		+ \vartheta\left(\mathbf{f}^{k+1}-\mathbf{f}^k\right)\right)
 * \f]
 * where
 * \f[
 * 	\mathbf{r}^k = \mathbf{f}^k - \mathrm{A}\mathbf{u}^k
 * \f]
 * and \f$ \mathbf{f}^{k+1} \f$ corresponds to the discretization of
 * \f$ f(t) \f$, where \f$ t \f$ is exactly the parameter passed as an input
 * to the method.
 */

/**
 * \fn void Problem::printSol(char const * filename)
 * \brief Prints in a file the array containing the solution
 *
 * This method recalls the function save_to_file of the Mesh class. In the
 * context of time-dependent problems, this function prints the solution
 * computed at the actual time.
 */

/**
 * \fn double Problem::evaluate(double const & t)
 * \brief Functional application on the solution array
 *
 * This will evaluate the functional FUN2 at the time t on the point of
 * the mesh and give back the values on the solution array u. In order to
 * make this method work, new versions of the functionals have been
 * written. Here is an example, but it is possible to find further
 * informations are in the related documentation.
 * \code
template<class EXACTSOL,class MESH,true>
class L2Error {
    private:
	template<class FUN>
	struct ActualTime {
	    double t;
	    FUN u;
	    ActualTime(double const & t): t(t) {}
	    double operator()(double const & x, double const & y) const {
		return u(x,y,t);
	    }
	};
    public:
        L2Error(MESH const * me, double const & t): me(me), uex(t) {}
	template<class T>
        double operator()(T const & unum, int n) const {
            return
                + me->prodL2(me->restriction(unum,n),
                             me->restriction(unum,n),n)

                + me->prodL2(me->restriction(uex,n),
                             me->restriction(uex,n),n)

                - 2 * me->prodL2(me->restriction(unum,n),
                                 me->restriction(uex,n),n);
        }
    private:
        MESH const * me;
        ActualTime<EXACTSOL> uex;
};
   \endcode
 */

/**
 * \fn void Problem::buildBC(std::valarray<double> & bc)
 * \brief Builds the boundary conditions on each edge
 *
 * The function supports either dirichlet and neumann boundary condition.
 * This method provide the discretization of the dirichlet and neumann
 * boundary terms into the array bc. In particular, for each point on the
 * border of the domain, it computes the imposed known value. If the border
 * is a dirichlet one, then the value of the solution is directly imposed
 * into the array bc. In the same way, also the value \f$ -\nabla\cdot u\f$
 * is discretized into the vector bc. The difference between the two
 * conditions will be exploited by the bcGen function, which will correctly
 * pach the known term vector b.
 */

/**
 * \fn double Problem::range(unsigned const & n) const
 * \brief Computes the range of the solution
 *
 * This method simply computes the range od the solution u at the actual
 * time. If the input parameter n is set to 0 the function returns the value
 * of the minimum. Instead, if the input parameter n is set to 1 (or to a
 * integer grater than 1) the function returns the maximum of the solution.
 * Owning a method like this is useful to allow the code genering by itself
 * visual representations of the problem's solution.
 */

/**
 * \fn void Problem::printPattern() const
 * \brief Prints a text file which contains the pattern of the stiffness matrix
 *
 * This method simply prints into a file called "matrixPattern.dat" the
 * structure of the stiffness matrix. It recalls the function, with the same
 * name, contained into the UMFMatrix class. In particular, it prints ones if
 * the considered value is different from zero. Then, using a visualization
 * software, it is possible to create an image of this pattern, for example
 * as with the Matlab command \c spy.
 */

/**
 * \fn void Problem::printAviScript(const char * const probname,
 *				    const char * const title,
 *				    int const & niter,
 *				    double const & min,
 *				    double const & max,
 *				    bool const & dim3) const
 * \brief Produces a bash script which generates a movie of the solution
 *
 * This method builds a bash script that works with GnuPlot and Mencoder.
 * It produces a video in avi format representing the solution at the
 * different temporal steps. The movie is obtained from the encoding of
 * a sequence of images produced with GnuPlot in a background mode. To
 * make this function work, it is necessary to have all the solution's
 * files named like \c probnameSol0000.dat, where probname is a label for
 * the specific problem and 0000 is a progressive number (of 4 digits)
 * which orders the several solutions. Moreover, all the solutions have
 * to be saved into a directory called \c probnameData.
 * The name of the output file is \c heat2mencoder.sh but, before executing
 * it with
 * \code
 * 	./heat2mencoder.sh
 * \endcode
 * it is necessary to give it the executable permissions through
 * \code
 * 	chmod 744 heat2mencoder.sh
 * \endcode
 *
 * \param probname A label for the kind of the problem
 * \param title The title which will show up in the header of the movie
 * \param niter The total number of solutions that have been collected
 * \param min Minimum value of the solution on the \f$[0,T]\f$ interval
 * \param max Maximum value of the solution on the \f$[0,T]\f$ interval
 * \param dim3 It if is true, a 3d movie is encoded. Otherwise the
 * 	  default value is false, in order to produce a 2d video.
 */

/**
 * \fn void Problem::printGnuplotScript(const char * const probname,
 * 					const char * const title,
 * 					double const & min,
 * 					double const & max,
 * 					bool const & dim3,
 * 					bool const & png) const
 * \brief Produces a bash script which generates pictures of the solution
 *
 * This method builds a bash script that works with GnuPlot. It produces
 * pictures of the problem's solution at different temporal steps. Those
 * pictures are obtained through a background call to GnuPlot and, to make
 * this function work, it is necessary that all the solution's file are
 * named like \c probnameSol0000.dat, where probname is a label for the
 * specific problem and 0000 is a progressive number (of 4 digits) which
 * orders the several solutions. Moreover, all the solution's files have
 * to be saved into a directory called \c probnameData. In this way, after
 * making the output file \c heat2gnuplot.sh executable through
 * \code
 * 	chmod 744 heat2gnuplot.sh
 * \endcode
 * it is possilbe to produce a single image calling
 * \code
 * 	./heat2gnuplot.sh NUM
 * \endcode
 * where \c NUM is the progressive number of the solution one wants to plot.
 * It's important to remark that the output pictures are in pdf, which is
 * an high quality vectorial format.
 *
 * \param probname A label for the kind of the problem
 * \param title The title which will show up in the header of the image
 * \param min Minimum value of the solution on the \f$[0,T]\f$ interval
 * \param max Maximum value of the solution on the \f$[0,T]\f$ interval
 * \param dim3 If it is true, a 3d graphic is plotted. Otherwise the
 *        default value is false, in order to produce a 2d figure
 * \param png If it is true the output will be in png format, otherwise
 *        (default) the code generates a pdf vector graphic.
 */

template<class BIL,class F>
class Problem {
    public:
        // Constructor
	Problem(Point<double> const & AP,
		Point<double> const & BP,
		int nx,
		int ny,
		int nnodes,
		double deltat,
		double theta,
		int nmax,
		double legDelta,
		double legEps,
		int order=1);

	// Initializes algebric data (without input)
	void initialize();

	// Initializes algebric data (with input)
	void initialize(UMFMatrix & A,
			UMFMatrix & M,
			UMFMatrix & K,
			std::valarray<double> & b,
			std::valarray<double> & bold,
			std::valarray<double> & u,
			std::valarray<double> & uold,
			std::valarray<double> & residual,
			std::valarray<double> & delta,
			std::valarray<double> const & bc);

        // Generates the stiffness matrix in the meaning of the vectors for UMF
        void matrixGen(UMFMatrix & A) const;

	// Generates the mass matrix in the meaning of the vectors for UMF
	void massMatrixGen(UMFMatrix & M) const;

	// Generates the vector of the initial data
	void initGen(std::valarray<double> & u) const;

        // Generates the vector of the known term
        void forceGen(std::valarray<double> & b) const;

	// Generates the vector of the known term when time dependent
	void forceGen(std::valarray<double> & b, double const & t) const;

	// Patches stiffness matrix to fit dirichlet boundary conditions
        void patchMatrix(UMFMatrix & A) const;

        // Applies boundary condition to the known term
        void bcGen(std::valarray<double> const & bc,
                   std::valarray<double> & b) const;

        // Changes the boudaries condition
	template<unsigned int border>
        void edgeBc(double const * const & newBcEdge,
                    std::valarray<double> & b) const;

	// Generates the K matrix for the time-advancing scheme
	void kappaMatrixGen(UMFMatrix & A, UMFMatrix & M, UMFMatrix & K);

	// Solves one step of the parabolic problem
	void solve();

	// Solves the parabolic problem with time-dependent forcing term
	void solve(double const & t);

	// Prints the solution into a file
	void printSol(char const * filename);

	// Builds the boundary conditions
	void buildBC(std::valarray<double> & bc);

	// Evaluates time dependent functional FUN2 at time t
	template<class FUN2>
	double evaluate(double const & t);

	// Computes the range of the solution
	double range(unsigned const & n) const;

	// Prints stiffness matrix pattern
	void printPattern() const;

	// Prints the script for Mencoder
	void printAviScript(const char * const probname,
			    const char * const title,
			    int const & niter,
			    double const & min,
			    double const & max,
			    bool const & dim3) const;

	// Prints the script for GnuPlot
	void printGnuplotScript(const char * const probname,
				const char * const title,
				double const & min,
				double const & max,
				bool const & dim3,
				bool const & png) const;

	// Mesh
	Mesh mesh;                     /**< \brief Mesh structure */
    private:
	// Algebric structures
	UMFMatrix _A;
	UMFMatrix _M;
	UMFMatrix _K;
	std::valarray<double> _b;
	std::valarray<double> _bold;
	std::valarray<double> _u;
	std::valarray<double> _uold;
	std::valarray<double> _residual;
	std::valarray<double> _delta;

	// Mesh parameter
	Point<double> _AP;              /**< \brief Bottom-right domain boundaries */
        Point<double> _BP;              /**< \brief Top-left domain boundaries */
        int _nx;                        /**< \brief Number of elements on the x-axis */
        int _ny;                        /**< \brief Number of elements on the y-axis */
	int _nnodes;                    /**< \brief Number of nodes in each element */
	double _deltat;			/**< \brief Length of the Time-step */
	double _theta;			/**< \brief Value for the relaxing parameter */
	int _sizeA;                     /**< \brief Problem dimension */
        int _nel;                       /**< \brief Total number of element */
        int _nNodesPerEl;               /**< \brief Number of nodes in an element */
        int _nBoundCond;                /**< \brief Number of nodes in the boundaries condition vector */

        // External problem functions
        BIL _bilinear;                  /**< \brief Bilinear form */

        // Functional
        F _functional;                  /**< \brief Functional */

        // Mesh polynomial
	mutable TensorProd _lagrange;	/**< \brief Legendre 1d polynomial */
        mutable TensorProd _lagrange1;  /**< \brief First legendre polynomial */
        mutable TensorProd _lagrange2;  /**< \brief Second legendre polynomial */

        // Boundaries conditions
	std::valarray<double> _bc;
        bctype _bt[4];                  /**< \brief Boundaries conditions type */
};

// ==========================================================
// End of Proper Header
// ==========================================================

template<class BIL,class F>
Problem<BIL,F>::Problem(Point<double> const & AP,
                        Point<double> const & BP,
                        int nx,
                        int ny,
                        int nnodes,
	    		double deltat,
			double theta,
                        int nmax,
                        double legDelta,
                        double legEps,
                        int order): mesh(nx,
                                         ny,
                                         nnodes,
                                         nmax,
                                         AP,
                                         BP,
                                         legDelta,
                                         legEps,
                                         order),
                                    _AP(AP),
                                    _BP(BP),
                                    _nx(nx),
                                    _ny(ny),
                                    _nnodes(nnodes),
				    _deltat(deltat),
				    _theta(theta),
                                    _bilinear(mesh),
                                    _functional(mesh) {
    // Computing boundaries conditions
    buildBC(_bc);

    // Mesh
    _nel = mesh.getnelem();
    _sizeA = mesh.getProblemDim();
    _nNodesPerEl = mesh.getnNodesPerEl();

    // Boundaries conditions
    _nBoundCond = mesh.getnbc();
    for(int i(0); i<4; ++i )  _bt[i] = _bilinear.bcond[i];
}

template<class BIL,class F>
void Problem<BIL,F>::initialize() {
    initialize(_A,_M,_K,_b,_bold,_u,_uold,_residual,_delta,_bc);
}

template<class BIL,class F>
void Problem<BIL,F>::initialize(UMFMatrix & A,
				UMFMatrix & M,
				UMFMatrix & K,
				std::valarray<double> & b,
				std::valarray<double> & bold,
				std::valarray<double> & u,
				std::valarray<double> & uold,
				std::valarray<double> & residual,
				std::valarray<double> & delta,
				std::valarray<double> const & bc) {
    A.resize_sizeA(_sizeA);
    M.resize_sizeA(_sizeA);
    K.resize_sizeA(_sizeA);
    b.resize(_sizeA);
    bold.resize(_sizeA);
    u.resize(_sizeA);
    uold.resize(_sizeA);
    residual.resize(_sizeA);
    delta.resize(_sizeA);
    matrixGen(A);
    massMatrixGen(M);
    forceGen(b,0);
    patchMatrix(A);
    bcGen(bc,b);
    initGen(u);
    kappaMatrixGen(A,M,K);
}

// Function for genering the matrix A
template<class BIL,class F>
void Problem<BIL,F>::matrixGen(UMFMatrix & A) const {
    int k;
    int I, J;

    // Computes matrix pattern
    std::list<int> * pattern;
    pattern = new std::list<int>[_sizeA];
    std::list<int>::const_iterator it;

    for(int idxCell(0); idxCell<_nel; ++idxCell) {
        for( int idxNod1(0); idxNod1<_nNodesPerEl; ++idxNod1 ) {
            for( int idxNod2(0); idxNod2<_nNodesPerEl; ++idxNod2 ) {
                pattern[mesh.getGlobalNum(idxCell,idxNod1)]
                       .push_front(mesh.getGlobalNum(idxCell,idxNod2));
            }
        }
    }

    // Fill the column indices array
    for( int i(0);i<_sizeA;++i ) {
        pattern[i].sort();
        pattern[i].unique();
    }

    A.getAp(0) = 0;
    for( int i(1); i<_sizeA; ++i ) {
        A.getAp(i) = A.getAp(i-1)+pattern[i-1].size();
    }

    int nNonZero = A.getAp(_sizeA)=A.getAp(_sizeA-1)+pattern[_sizeA-1].size();
    A.resize_sizeNZ(nNonZero);

    int idx(0);
    for(int i(0);i<_sizeA;++i) {
        idx=0;
        for(it=pattern[i].begin();it!=pattern[i].end();++it) {
            A.getAi(idx+A.getAp(i))=*it;
            ++idx;
        }
    }
    delete[] pattern;

    for(int idxCell(0); idxCell<_nel; ++idxCell )  {
        for(int idxNod1(0); idxNod1<_nNodesPerEl; ++idxNod1 )  {
            _lagrange1 = mesh.lagrange( idxNod1%_nnodes, idxNod1/_nnodes );
            I = mesh.getGlobalNum(idxCell,idxNod1);

            for(int idxNod2(0); idxNod2<_nNodesPerEl; ++idxNod2)  {
                _lagrange2 = mesh.lagrange( idxNod2%_nnodes,idxNod2/_nnodes );

                J = mesh.getGlobalNum(idxCell,idxNod2);
                k = A.getAp(J);
                while( A.getAi(k)!= I ) {
                    ++k;
                }
                assert( k<A.getAp(J+1) );

                // J is the row index here:
                A.getAx(k) += _bilinear(_lagrange1,_lagrange2,idxCell);
            }
        }
    }
}

// Function for genering the matrix M
template<class BIL,class F>
void Problem<BIL,F>::massMatrixGen(UMFMatrix & M) const {
    int k;
    int I, J;

    // Computes matrix pattern
    std::list<int> * pattern;
    pattern = new std::list<int>[_sizeA];
    std::list<int>::const_iterator it;

    for(int idxCell(0); idxCell<_nel; ++idxCell) {
        for( int idxNod1(0); idxNod1<_nNodesPerEl; ++idxNod1 ) {
            for( int idxNod2(0); idxNod2<_nNodesPerEl; ++idxNod2 ) {
                pattern[mesh.getGlobalNum(idxCell,idxNod1)]
                       .push_front(mesh.getGlobalNum(idxCell,idxNod2));
            }
        }
    }

    // Fill the column indices array
    for( int i(0);i<_sizeA;++i ) {
        pattern[i].sort();
        pattern[i].unique();
    }

    M.getAp(0) = 0;
    for( int i(1); i<_sizeA; ++i ) {
        M.getAp(i) = M.getAp(i-1)+pattern[i-1].size();
    }

    int nNonZero = M.getAp(_sizeA)=M.getAp(_sizeA-1)+pattern[_sizeA-1].size();
    M.resize_sizeNZ(nNonZero);

    int idx(0);
    for(int i(0);i<_sizeA;++i) {
        idx=0;
        for(it=pattern[i].begin();it!=pattern[i].end();++it) {
            M.getAi(idx+M.getAp(i))=*it;
            ++idx;
        }
    }
    delete[] pattern;

    for(int idxCell(0); idxCell<_nel; ++idxCell )  {
        for(int idxNod1(0); idxNod1<_nNodesPerEl; ++idxNod1 )  {
            _lagrange1 = mesh.lagrange( idxNod1%_nnodes, idxNod1/_nnodes );
            I = mesh.getGlobalNum(idxCell,idxNod1);

            for(int idxNod2(0); idxNod2<_nNodesPerEl; ++idxNod2)  {
                _lagrange2 = mesh.lagrange( idxNod2%_nnodes,idxNod2/_nnodes );

                J = mesh.getGlobalNum(idxCell,idxNod2);
                k = M.getAp(J);
                while( M.getAi(k)!= I ) {
                    ++k;
                }
                assert( k<M.getAp(J+1) );

                // J is the row index here:
                M.getAx(k) += mesh.prodL2(_lagrange1,_lagrange2,idxCell);
            }
        }
    }
}

// Function for genering initial data
template<class BIL,class F>
void Problem<BIL,F>::initGen(std::valarray<double> & u) const {
    typedef typename BIL::initial initData;
    initData u0;
    mesh.interpolate(u0, u);
}

// Function for genering the forcing term
template<class BIL,class F>
void Problem<BIL,F>::forceGen(std::valarray<double> & b) const  {
    int I;

    for(int idxCell(0); idxCell<_nel; ++idxCell)  {
        for(int idxNod1(0); idxNod1<_nNodesPerEl; ++idxNod1 )  {

            _lagrange1 = mesh.lagrange( idxNod1%_nnodes, idxNod1/_nnodes );
            I = mesh.getGlobalNum(idxCell,idxNod1);
            assert(I>=0 && I<_sizeA);

            b[I] += _functional(_lagrange1,idxCell);
        }
    }
}

// Function for genering the forcing term when time dependent
template<class BIL,class F>
void Problem<BIL,F>::forceGen(std::valarray<double> & b, double const & t) const  {
    F functional(mesh, t);
    int I;

    for(int idxCell(0); idxCell<_nel; ++idxCell)  {
        for(int idxNod1(0); idxNod1<_nNodesPerEl; ++idxNod1 )  {

            _lagrange1 = mesh.lagrange( idxNod1%_nnodes, idxNod1/_nnodes );
            I = mesh.getGlobalNum(idxCell,idxNod1);
            assert(I>=0 && I<_sizeA);

            b[I] += functional(_lagrange1,idxCell);
        }
    }
}

// Boundary condition function generator
template<class BIL,class F>
void Problem<BIL,F>::patchMatrix(UMFMatrix & A) const  {
    int I;

    std::vector<int> elements;
    std::vector<int> nodes;

    std::vector<int>::const_iterator element;
    std::vector<int>::const_iterator node;

    for(int border=0;border<4;++border) {
        if(_bt[border]==dirichlet) {
            mesh.getBordElements(border,elements);
            for(element=elements.begin();element<elements.end();++element) {
                mesh.getGlobBordNodes(border,*element,nodes);
                for(node=nodes.begin();node<nodes.end();++node) {
                    I=*node;
                    assert(I>=0 && I<_sizeA);
                    for(int row=A.getAp(I);row<A.getAp(I+1);++row) {
                        A.getAx(row)=( (A.getAi(row)==I) ? 1.0 : 0.0 );
                    }
                }
            }
        }
    }
}

template<class BIL,class F>
void Problem<BIL,F>::bcGen(std::valarray<double> const & bc,
			   std::valarray<double> & b) const {
    int k(0);

    edgeBc<0>(&bc[k],b);
    k+=mesh.getnNodesPerBord(0);

    edgeBc<1>(&bc[k],b);
    k+=mesh.getnNodesPerBord(1);

    edgeBc<2>(&bc[k],b);
    k+=mesh.getnNodesPerBord(2);

    edgeBc<3>(&bc[k],b);
    k+=mesh.getnNodesPerBord(3);
}

template<class BIL,class F>
template<unsigned int border>
void Problem<BIL,F>::edgeBc(double const * const & newBcEdge,
                            std::valarray<double> & b) const {
    int I;

    std::vector<int> elements;
    std::vector<int> nodes;
    std::vector<int>::const_iterator element;
    std::vector<int>::const_iterator node;

    int k=0;

    if (_bt[border]==dirichlet) {
	mesh.getBordElements(border,elements);
	for(element=elements.begin();element<elements.end();++element) {
	    mesh.getGlobBordNodes(border,*element,nodes);
	    for(node=nodes.begin();node<nodes.end();++node) {
		I=*node;
		assert(I>=0 && I<_sizeA);
		b[I]=newBcEdge[mesh.getBordNum(border,I)];
	    }
	}
    } else if (_bt[border]==neumann) {
	mesh.getBordElements(border,elements);
	for(element=elements.begin();element<elements.end();++element) {
	    mesh.getGlobBordNodes(border,*element,nodes);
	    for(node=nodes.begin();node<nodes.end();++node) {
		_lagrange=mesh.lagrange(k);
		I=*node;
		assert(I>=0 && I<_sizeA);
		typedef typename BIL::template on<border> OnNeuBorder;
		OnNeuBorder _bcneu;
		b[I]+=mesh.prodL2border(mesh.restriction(_bcneu,*element),
			_lagrange,*element,border);
		++k;
	    }
	    k=0;
	}
    }
}

template<class BIL,class F>
void Problem<BIL,F>::kappaMatrixGen(UMFMatrix & A,
				    UMFMatrix & M,
				    UMFMatrix & K) {
    K = _theta * A + (1./_deltat) * M;
    K.factorize();
}

template<class BIL,class F>
void Problem<BIL,F>::solve() {
    _uold = _u;
    _A.multiply(_uold, _residual);
    _residual = _b - _residual;
    _K.solve(_residual, _delta);
    _u = _uold + _delta;
}

template<class BIL,class F>
void Problem<BIL,F>::solve(double const & t) {
    _uold = _u;
    _bold = _b;
    forceGen(_b, t);
    bcGen(_bc, _b);
    _A.multiply(_uold, _residual);
    _residual = _theta*_b + (1-_theta)*_bold - _residual;
    _K.solve(_residual, _delta);
    _u = _uold + _delta;
}

template<class BIL,class F>
void Problem<BIL,F>::printSol(char const * filename) {
    mesh.save_to_file(_u, filename);
}

template<class BIL,class F>
template<class FUN2>
double Problem<BIL,F>::evaluate(double const & t) {
    FUN2 J(&mesh, t);
    double sum=0.e0;
    for(int n=0;n<_nel;++n) {
	sum+=J(_u,n);
    }
    return sum;
}

template<class BIL,class F>
void Problem<BIL,F>::buildBC(std::valarray<double> & bc) {
    Point<double> P;

    // Number of nodes
    int _nnx = mesh.getnx()*(mesh.getnnodes()-1) + 1;
    int _nny = mesh.getny()*(mesh.getnnodes()-1) + 1;

    bc.resize(2*_nnx + 2*_nny);
    int k(0);

    std::vector<int> bnodes;
    std::vector<int>::const_iterator it;

    // Edge 0
    mesh.getDomainBordNodes(0, bnodes);
    for(it = bnodes.begin(); it<bnodes.end(); ++it) {
	P = mesh.getNodePosition(*it);
	typedef typename BIL::template on<0> OnLowerEdge;
	OnLowerEdge _down;
	bc[k++] = _down(P[0], P[1]);
    }

    // Edge 1
    mesh.getDomainBordNodes(1, bnodes);
    for(it = bnodes.begin(); it<bnodes.end(); ++it) {
	P = mesh.getNodePosition(*it);
	typedef typename BIL::template on<1> OnRightEdge;
	OnRightEdge _right;
	bc[k++] = _right(P[0], P[1]);
    }

    // Edge 2
    mesh.getDomainBordNodes(2, bnodes);
    for(it = bnodes.begin(); it<bnodes.end(); ++it) {
	P = mesh.getNodePosition(*it);
	typedef typename BIL::template on<2> OnUpperEdge;
	OnUpperEdge _up;
	bc[k++] = _up(P[0], P[1]);
    }

    // Edge 3
    mesh.getDomainBordNodes(3, bnodes);
    for(it = bnodes.begin(); it<bnodes.end(); ++it) {
	P = mesh.getNodePosition(*it);
	typedef typename BIL::template on<3> OnLeftEdge;
	OnLeftEdge _left;
	bc[k++] = _left(P[0], P[1]);
    }

    // Postcondition
    assert(k == 2*_nnx + 2*_nny);
}

template<class BIL,class F>
double Problem<BIL,F>::range(unsigned const & n) const {
    return ( (n==0) ? _u.min() : _u.max() );
}

template<class BIL,class F>
void Problem<BIL,F>::printPattern() const {
    _A.printPattern("matrixPattern.dat");
}

template<class BIL,class F>
void Problem<BIL,F>::printAviScript(const char * const probname,
				    const char * const title,
				    int const & niter,
				    double const & min,
				    double const & max,
				    bool const & dim3=false) const {
    int horpixel(0);
    int verpixel(0);
    double xlength = _BP.get(0) - _AP.get(0);
    double ylength = _BP.get(1) - _AP.get(1);
    double ratio = xlength / ylength;
    
    if (ratio >= 1) {
	horpixel = 800;
	verpixel = int(800 / ratio);
    } else {
	verpixel = 800;
	horpixel = int(800 * ratio);
    }

    std::ofstream encoder;
    encoder.open("heat2mencoder.sh", std::ios_base::out);

    encoder << "#!/bin/bash" << std::endl << std::endl
	    << "name=" << probname << std::endl
	    << "dir=" << probname << "Data" << std::endl
	    << "cd $dir" << std::endl
	    << "for n in {0.." << niter << "}" << std::endl
	    << "do" << std::endl
	    << "printf -v index '%04d' $n" << std::endl
	    << "gnuplot << EOF" << std::endl
	    << "set terminal jpeg medium size " << horpixel << "," << verpixel << std::endl
	    << "set output '" << probname << "Sol$index.jpg'" << std::endl
	    << "set title '" << title << "'" << std::endl
	    << "set xrange [" << _AP.get(0) << ":" << _BP.get(0) << "]" << std::endl
	    << "set yrange [" << _AP.get(1) << ":" << _BP.get(1) << "]" << std::endl
	    << "set zrange [" << min << ":" << max << "]" << std::endl
	    << "set key off" << std::endl;
    if (!dim3) encoder << "set pm3d map" << std::endl;
    else encoder << "set pm3d at s" << std::endl
		 << "unset surf" << std::endl;
    encoder << "set palette rgbformulae 21,22,23" << std::endl
	    << "set cbrange [" << min << ":" << max << "]" << std::endl
	    << "time = $index * " << _deltat << std::endl;
    if (!dim3) encoder << "set label 1 at " << _BP.get(0)-1 << ", " << _AP.get(1)-0.5 << std::endl;
    else encoder << "set label 1 at 5, -9" << std::endl;
    encoder << "set label 1 sprintf('Time: %.2f sec', time)" << std::endl
	    << "splot '" << probname << "Sol$index.dat" << std::endl
	    << "EOF" << std::endl
	    << "done" << std::endl
	    << "mencoder 'mf://*.jpg' -mf fps=25 -o $name.avi -ovc lavc -lavcopts vcodec=mpeg4:vbitrate=800" << std::endl
	    << "rm -rf *.jpg" << std::endl
	    << "mv $name.avi ../";

    encoder.close();
}

template<class BIL,class F>
void Problem<BIL,F>::printGnuplotScript(const char * const probname,
					const char * const title,
					double const & min,
					double const & max,
					bool const & dim3=false,
					bool const & png=false) const {
    int horpixel(0), verpixel(0);
    double horinch(0), verinch(0);
    double xlength = _BP.get(0) - _AP.get(0);
    double ylength = _BP.get(1) - _AP.get(1);
    double ratio = xlength / ylength;
    
    if (ratio >= 1) {
	horinch = 5.;
	horpixel = 800;
	verinch = 5. / ratio;
	verpixel = int(800 / ratio);
    } else {
	verinch = 5.;
	verpixel = 800;
	horinch = 5. * ratio;
	horpixel = int(800 * ratio);
    }

    std::ofstream gnuplot;
    gnuplot.open("heat2gnuplot.sh", std::ios_base::out);

    gnuplot << "#!/bin/bash" << std::endl << std::endl
	    << "name=" << probname << std::endl
	    << "dir=" << probname << "Data" << std::endl
	    << "cd $dir" << std::endl
	    << "printf -v index '%04d' $1" << std::endl
	    << "gnuplot << EOF" << std::endl;
    if (!png) gnuplot << "set terminal pdf enhanced color fsize 8 size " << horinch << "," << verinch << std::endl
	<< "set output '" << probname << "Sol$index.pdf'" << std::endl;

    else gnuplot << "set terminal png medium size " << horpixel << "," << verpixel << std::endl
	<< "set output '" << probname << "Sol$index.png'" << std::endl;
    gnuplot << "set title '" << title << "'" << std::endl
	    << "set xrange [" << _AP.get(0) << ":" << _BP.get(0) << "]" << std::endl
	    << "set yrange [" << _AP.get(1) << ":" << _BP.get(1) << "]" << std::endl
	    << "set zrange [" << min << ":" << max << "]" << std::endl
	    << "set key off" << std::endl;
    if (!dim3) gnuplot << "set pm3d map" << std::endl;
    else gnuplot << "set pm3d at s" << std::endl
		 << "unset surf" << std::endl;
    gnuplot << "set palette rgbformulae 21,22,23" << std::endl
	    << "set cbrange [" << min << ":" << max << "]" << std::endl
	    << "splot '" << probname << "Sol$index.dat'" << std::endl
	    << "EOF" << std::endl;
    if (!png) gnuplot << "mv *.pdf ../";
    else gnuplot << "mv *.png ../";
    gnuplot.close();
}

} // end namespace

#endif
